home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / Lyapunov / LyapunovView.m < prev    next >
Encoding:
Text File  |  1992-08-01  |  14.9 KB  |  502 lines
  1. // Note: you can't yet render only parts of the fractal.
  2. // Eventually the engine needs to be moved to a separate
  3. // object from this view so that several engines -- each
  4. // with it's own thread -- can be started up and all feed
  5. // into this view object.
  6.  
  7. // Adapted from lyap.c by Don Yacktman
  8. // This object embodies the lyapunov calculation code.
  9. // Use the methods provided to set the various options.
  10. // You can take an already-calculated fractal and apply
  11. // a color map to it through methods provided.  Just pass
  12. // it a color map object that you have created.
  13.  
  14. // Here's the copyright from the original UNIX program "lyap",
  15. // from which I derived much of the engine.
  16. // The original program and an executable and output viewer for
  17. // the NeXT are provided in a separate directory in this distribution.
  18. // As the notice states, you cannot redistribute this program, or
  19. // anything derived from it without including all sources, including
  20. // your own changes!  
  21.  
  22. /*
  23.  * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
  24.  * This program is Copyright (C) 1991, Garrett A. Wollman.  This
  25.  * program may be modified and distributed for any purpose and without
  26.  * fee, provided that this notice remains on all such copies
  27.  * unaltered.  Binary distributions not including this source file are
  28.  * prohibited.  Modified versions must be marked with the name of the
  29.  * modifier and the date of modification.
  30.  *
  31.  * Under no circumstances shall Garrett A. Wollman or the University
  32.  * of Vermont and State Agricultural College be held liable for any
  33.  * damages resulting from the use or misuse of this program, whether
  34.  * the author is aware of such possibility or not.  This program is
  35.  * warranted solely to occupy disk space.
  36.  *
  37.  * I'm sorry for all this legal garbage, but it is sadly necessary
  38.  * these days.
  39.  */
  40.  
  41. /* Generated by Interface Builder */
  42.  
  43. #import "LyapunovView.h"
  44. #import "ColorMapView.h"
  45. #import "ColorMap.h"
  46. #import <libc.h>
  47. #import <stdlib.h>
  48. #import <stdio.h>
  49. #import <strings.h>
  50. #import <appkit/tiff.h>
  51. #import <appkit/color.h>
  52. #import <appkit/graphics.h>
  53. #import <appkit/Control.h>
  54. #import <appkit/SavePanel.h>
  55. #import <appkit/OpenPanel.h>
  56. #import <appkit/Matrix.h>
  57. #import <appkit/NXBitmapImageRep.h>
  58. #import <appkit/NXColorWell.h>
  59. #import <dpsclient/psops.h>
  60. #import <appkit/Application.h>
  61. #import <dpsclient/wraps.h>
  62. #import <sys/file.h>
  63. #import <streams/streams.h>
  64. #import <sys/types.h>
  65. #import <sys/uio.h>
  66.  
  67. #define FABS(k) (((k) < 0) ? (-(k)) : (k))    // faster than 2.1 trap
  68.  
  69. double log(x)    // to replace trap in 2.x OS that slows down '040 machines
  70. double x;        // taken verbatim from libjv.a
  71. {
  72. /*
  73.     Routine for log(x) in IEEE double precision.
  74.     Note that we have to hack around in the exponent!
  75.     Steps: 1 x = y*2^n, with y between sqrt(2) and 1/sqrt(2)
  76.            2 convert to z = (y-1)/(y+1)
  77.            3 ln = 2*z*sum(i=0,9, z^(2*i)/(2*i+1)) + n*ln(2)
  78.     Routine made by J.A.M. Vermaseren 20-apr-1992.
  79. */
  80.     double scr, y, z, z2, sum, sqrt();
  81.     short n;
  82.     if ( x < 0 ) {
  83.         printf("log of a negative number: %d\n",x);
  84.         return (-1);
  85.     }
  86.     scr = x;
  87.     n = ( (*((short *)(&scr))) & 0x7FF0 ) - 0x3FF0;
  88.     *((short *)(&scr)) -= n;
  89.     n >>= 4;
  90.     if ( scr > 1.4142135623730950488 ) {
  91.         *((short *)(&scr)) -= 0x10;
  92.         n++;
  93.     }
  94.     y = scr;
  95.     z = (y-1.)/(y+1.);
  96.     z2 = z*z;
  97.     sum = 2.*z*(0.99999999999999999886+z2*(
  98.                 0.33333333333333825805+z2*(
  99.                 0.19999999999649522967+z2*(
  100.                 0.14285714380662481221+z2*(
  101.                 0.11111098494595009814+z2*(
  102.                 0.09091817553294820737+z2*(
  103.                 0.07656220982777783836+z2*(
  104.                 0.07405280893003482131))))))))
  105.          +n*0.6931471805599453094;
  106.     return(sum);
  107. }
  108.  
  109. @implementation LyapunovView
  110.  
  111. - initFrame:(const NXRect *)frm        // designated initializer for a view
  112. {
  113.     [super initFrame:frm];
  114.     xPos = 0.0; yPos = 0.0;
  115.     xScale = 4.0; yScale = 4.0;
  116.     initial = 0.5; minExp = -5.0;
  117.     dwell = 2000;  settle = 600;
  118.     firstInit = YES; positive = NO;
  119.     return self;
  120. }
  121.  
  122. - appDidInit:sender                // final initialization
  123. {            // set up default parameters on control panel
  124.     activeMap = [activeMapView colorMap];
  125.     
  126.     // set up default position
  127.     [[posMatrix findCellWithTag:0] setFloatValue:xPos];
  128.     [[posMatrix findCellWithTag:1] setFloatValue:yPos];
  129.     
  130.     // set up default scale
  131.     [[scaleMatrix findCellWithTag:0] setFloatValue:xScale];
  132.     [[scaleMatrix findCellWithTag:1] setFloatValue:yScale];
  133.     
  134.     // set up upper r.h. corner
  135.     [[sizeMatrix findCellWithTag:0] setFloatValue:(xPos + xScale)];
  136.     [[sizeMatrix findCellWithTag:1] setFloatValue:(yPos + yScale)];
  137.     [[windowSizeMatrix findCellWithTag:0] setFloatValue:(XSIZE)];
  138.     [[windowSizeMatrix findCellWithTag:1] setFloatValue:(YSIZE)];
  139.     
  140.     // set up default initial value
  141.     [[boundsMatrix findCellWithTag:1] setFloatValue:initial];
  142.     [[boundsMatrix findCellWithTag:0] setFloatValue:minExp];
  143.     
  144.     // set up default depth
  145.     [[depthMatrix findCellWithTag:0] setIntValue:dwell];
  146.     [[depthMatrix findCellWithTag:1] setIntValue:settle];
  147.     
  148.     // set up polarity
  149.     [polarityMatrix selectCellWithTag:positive];
  150.  
  151.     // set up default pattern
  152.     if (firstInit) [[pattern findCellWithTag:0] setStringValue:"aaab"];
  153.     firstInit = NO;
  154.     [saveMenuButton setEnabled:NO];
  155.     [reApplyButton setEnabled:NO];
  156.     return self;
  157. }
  158.  
  159.  
  160.  
  161. // Lyapunov exponent and coloring calculation
  162. // original function by Ed Kubaitis, used by permission
  163.  
  164. // Speedup...
  165. // Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
  166. // is, by a principle from high-school algebra, the same as
  167. // $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
  168. // by unrolling the dwell loop by four (hence the rounding in main())
  169. // and accumulating by multiplication before taking the log.  Since
  170. // log() is extremely expensive, this should lead to a considerable
  171. // speedup, while still allowing complete flexibility in selecting
  172. // dwell values.
  173. // Presumably, one might unroll this loop even further, with a 
  174. // concomitant increase in textual complexity.
  175.  
  176. - (int)lya:(double)a :(double)b
  177. {
  178.     double t = 0;
  179.     double r, lya;
  180.     int d;
  181.     double acc;
  182.     double x;
  183.  
  184.     x = initial;        /* initialize x */
  185.  
  186.     for(d=0; d < settle; d++) {
  187.         r = (patternNum[d]) ? b : a;
  188.         x = r*x*(1-x);
  189.     }
  190.  
  191.     for(d=0; d < dwell; d+= 4) {
  192.         r = (patternNum[d]) ? b : a;
  193.         x = r * x * (1-x);
  194.         acc = r - 2*r*x;
  195.  
  196.         r = (patternNum[d+1]) ? b : a;
  197.         x = r * x * (1-x);
  198.         acc *= r - 2*r*x;
  199.  
  200.         r = (patternNum[d+2]) ? b : a;
  201.         x = r * x * (1-x);
  202.         acc *= r - 2*r*x;
  203.  
  204.         r = (patternNum[d+3]) ? b : a;
  205.         x = r * x * (1-x);
  206.         t += log(FABS(acc * (r - 2*r*x)));
  207.  
  208.         if(FABS(x-0.5) < 1e-20) {
  209.             d += 3;
  210.             break;
  211.         }
  212.     }
  213.  
  214.     lya = (t * 1.442695041)/d;
  215.  
  216.     if (positive) lya *= -1;
  217.     if (lya < minExp) lya = minExp;
  218.  
  219.     if (lya < 0) {
  220.         return (int)(lya / minExp * (double)NUMCOLORS);
  221.     } else {
  222.         return 0;
  223.     }
  224. }
  225.  
  226. - go:sender                // begin calculation of image
  227. {
  228.     register int pix, y, x;
  229.     double aa, bb;
  230.     int r, g, b;
  231.     unsigned char c;
  232.     NXRect pixel, line;
  233.     
  234.     pixel.size.width = 1.0;
  235.     pixel.size.height = 1.0;
  236.     line.size.height = 1.0;
  237.     line.size.width = XSIZE;
  238.     line.origin.x = 0.0;
  239.     
  240.     [self newParam:self];    // be sure we're up to date
  241.     
  242.     // main calc loop: calc line & draw it.
  243.     [self lockFocus];
  244.     for (y=0; y<YSIZE; y++) {
  245.         bb = yPos + y*yScale/YSIZE; //  get physical y-coord
  246.         for (x=0; x<XSIZE; x++) {
  247.             aa = xPos + x*xScale/XSIZE; //  get physical x-coord
  248.             space[x][y] = [self lya:aa :bb]; // Lyapunov exponent approx.
  249.             // put into pixel map w/color        
  250.             pix = (x+(YSIZE-y-1)*XSIZE)*RGB;
  251.             c = (unsigned char)space[x][y];
  252.             [activeMap colorFor:c :&r :&g :&b];
  253.             pixels[pix] = r;  // pix+RED = pix because RED=0
  254.             pixels[pix+GREEN] = g;
  255.             pixels[pix+BLUE] = b;
  256.         }
  257.         // plot a line of it
  258.         line.origin.y = y;
  259.         NXImageBitmap(&line, XSIZE, 1, 8, 3,
  260.             NX_MESHED, NX_COLORMASK,
  261.             &pixels[(int)((YSIZE-y-1)*XSIZE*RGB+0.5)],
  262.         NULL, NULL, NULL, NULL);
  263.         [window flushWindow];
  264.         NXPing();
  265.     }
  266.     [self unlockFocus];
  267.     [saveMenuButton setEnabled:YES];
  268.     [reApplyButton setEnabled:YES];
  269.     return self;
  270. }
  271.  
  272. - reApply:sender        // re do coloring on plot
  273. {
  274.     register int x, y, pix;
  275.     register unsigned char c;
  276.     int r, g, b;
  277.     NXRect pixel, line;
  278.     
  279.     [self newParam:self];    // be sure we're up to date
  280.     // set up plotting bounds
  281.     pixel.size.width = 1; pixel.size.height = 1;
  282.     line.size.height = 1; line.size.width = XSIZE;
  283.     line.origin.x = 0;
  284.     // re-apply colors
  285.     [self lockFocus];
  286.     for (y=0; y<YSIZE; y++) {
  287.         for (x=0; x<XSIZE; x++) {
  288.             pix = (x+(YSIZE-y-1)*XSIZE)*RGB;
  289.             c = (unsigned char)space[x][y];
  290.             [activeMap colorFor:c :&r :&g :&b];
  291.             pixels[pix] = r;  // pix+RED = pix because RED=0
  292.             pixels[pix+GREEN] = g;
  293.             pixels[pix+BLUE] = b;
  294.         }
  295.         // plot a line of it
  296.         line.origin.y = y;
  297.         NXImageBitmap(&line, XSIZE, 1, 8, 3,
  298.             NX_MESHED, NX_COLORMASK,
  299.             &pixels[(int)((YSIZE-y-1)*XSIZE*RGB+0.5)],
  300.             NULL, NULL, NULL, NULL);
  301.         [window flushWindow];
  302.         NXPing();
  303.     } 
  304.     [self unlockFocus];   
  305.     [saveMenuButton setEnabled:YES];
  306.     [reApplyButton setEnabled:YES];
  307.     return self;
  308. }
  309.  
  310. - newParam:sender            // set up new parameters for calc
  311. {
  312.     int count, len;
  313.     
  314.     // get x-y coordinates
  315.     xPos = [[posMatrix findCellWithTag:0] floatValue];
  316.     yPos = [[posMatrix findCellWithTag:1] floatValue];
  317.     
  318.     // get x-y scaling
  319.     xScale = [[scaleMatrix findCellWithTag:0] floatValue];
  320.     yScale = [[scaleMatrix findCellWithTag:1] floatValue];
  321.     
  322.     // get initial value for each point
  323.     minExp  = [[boundsMatrix findCellWithTag:0] floatValue];
  324.     initial = [[boundsMatrix findCellWithTag:1] floatValue];
  325.     
  326.     // get depth
  327.     settle = [[depthMatrix findCellWithTag:1] intValue];
  328.     dwell = [[depthMatrix findCellWithTag:0] intValue];
  329.     if (dwell % UNRAVEL) dwell += UNRAVEL - (dwell % UNRAVEL); // round off
  330.     if (dwell < settle) {    // sanity check
  331.         if (settle % UNRAVEL) settle += UNRAVEL - (settle % UNRAVEL);
  332.         dwell = settle;
  333.     }
  334.     
  335.     // get pattern
  336.     patternString = [[pattern findCellWithTag:0] stringValue];
  337.     len = strlen(patternString);
  338.     if (patternNum) free(patternNum);
  339.     patternNum = (int *)malloc((dwell + 1) * sizeof (int));
  340.     if(!patternNum) {
  341.         fprintf(stderr, "Cannot allocate space for Markus vector!\n");
  342.     }
  343.     for(count=0; count <= dwell; count++) {
  344.         patternNum[count] = ((patternString[count % len] == 'a') ||
  345.             (patternString[count % len] == 'x')) ? 0 : 1;
  346.     }
  347.     patternLength = count;
  348.     
  349.     positive = [[polarityMatrix selectedCell] tag];
  350.     
  351.     [self appDidInit:self];
  352.     return self;
  353. }
  354.  
  355. - drawSelf:(NXRect *)rects :(int)rectCount    // redraws the screen.
  356. {
  357.     PSsetgray(NX_DKGRAY);
  358.     NXRectFill(&bounds);
  359.     return self;
  360. }
  361.  
  362. - mouseDown:(NXEvent *)e    // stolen from MandelView.m
  363. /*
  364.  * We implement the mouseDown method so the user can sweep out a section of
  365.  * the view and select that as the current x,y,scale coordinates.
  366.  */
  367. {
  368.   int looping = YES;
  369.   int oldMask;
  370.   NXRect bbox;
  371.   NXPoint startPoint, currPoint;
  372.   double newX, newY, newDX, newDY;
  373.   NXCoord aspectRatio;
  374.   float mvX = xPos;
  375.   float mvY = yPos;
  376.   float mvDX = xScale;
  377.   float mvDY = yScale;
  378.   
  379.   aspectRatio = bounds.size.height / bounds.size.width;
  380.  
  381.   oldMask =  [window addToEventMask:NX_MOUSEDRAGGEDMASK];
  382.   startPoint = e->location;
  383.   [self convertPoint:&startPoint fromView:nil];
  384.   NXSetRect(&bbox,startPoint.x,startPoint.y,0.0,0.0);
  385.   [self lockFocus];
  386.   while (looping) {
  387.     e=[NXApp getNextEvent:NX_MOUSEUPMASK | NX_MOUSEDRAGGEDMASK];
  388.     currPoint = e->location;
  389.     [self convertPoint:&currPoint fromView:nil];
  390.     bbox.size.width = 2*(currPoint.x - startPoint.x);
  391.     bbox.size.height = 2*(currPoint.y - startPoint.y);
  392.     /* Normalize bbox to always have positive width and height */
  393.     if (bbox.size.width < 0) bbox.size.width = -bbox.size.width;
  394.     if (bbox.size.height < 0) bbox.size.height = -bbox.size.height;
  395.     /*
  396.      * constrain the box to have the aspect ratio of the view.  Choose
  397.      * whichever dimension is closer to the desired result.
  398.      */
  399.     if ((bbox.size.height/bbox.size.width) > aspectRatio) {
  400.       bbox.size.height = bbox.size.width * aspectRatio;
  401.     } else {
  402.       bbox.size.width = bbox.size.height / aspectRatio;
  403.     }
  404.     /* The startPoint is always at the center of the bbox */
  405.     bbox.origin.x = startPoint.x - .5 * bbox.size.width;
  406.     bbox.origin.y = startPoint.y - .5 * bbox.size.height;
  407.     PSnewinstance();
  408.     if (looping = (e->type == NX_MOUSEDRAGGED)) {
  409.       PSsetinstance(YES);
  410.       PSsetgray(NX_WHITE);
  411.       NXFrameRect(&bbox);
  412.       PSsetinstance(NO);
  413.     }
  414.   }
  415.   [self unlockFocus];
  416.   [window setEventMask:oldMask];
  417.   if ((bbox.size.width > 0) && (bbox.size.height > 0)) {
  418.     /*
  419.      * At this point, bbox is in window coordinates.  Set new controller
  420.      * parameters based on this view's coordinates and the bounding box.
  421.      */
  422.     newDX = (bbox.size.width*mvDX)/bounds.size.width;
  423.     newDY = (bbox.size.height*mvDY)/bounds.size.height;
  424.     newX = ((startPoint.x + bounds.origin.x) / bounds.size.width);
  425.     newX = newX * mvDX + mvX - newDX/2;
  426.     newY = ((startPoint.y + bounds.origin.y) / bounds.size.height);
  427.     newY = newY * mvDY + mvY - newDY/2;
  428.  
  429.     /*
  430.      * Note that we only update the text fields -- we
  431.      * don't update the internal instance variables.  If we were to update
  432.      * the internal ivs, then the user would only get one chance with the
  433.      * mouse down method because subsequent mouse-downs would work in terms
  434.      * of the new coordinate system.
  435.      */
  436.     [[posMatrix findCellWithTag:0] setFloatValue:newX];
  437.     [[posMatrix findCellWithTag:1] setFloatValue:newY];
  438.     [[scaleMatrix findCellWithTag:0] setFloatValue:newDX];
  439.     [[scaleMatrix findCellWithTag:1] setFloatValue:newDY];
  440.   }
  441.   return self;
  442. }
  443.  
  444. - saveAsTIFF:sender;            // save image in a .tiff file
  445. {
  446.     if ([[SavePanel new] runModalForDirectory:[[OpenPanel new] directory]
  447.             file:"UNTITLED.tiff"]) {
  448.         [self saveTIFF:[[SavePanel new] filename]];
  449.     }
  450.     return self;
  451. }
  452.  
  453. - (BOOL)saveTIFF:(const char *)name
  454. {
  455.     id imageRep;
  456.     int fd;
  457.     NXStream *stream;
  458.  
  459.     imageRep = [[NXBitmapImageRep alloc] initData:pixels
  460.         pixelsWide:(int)(XSIZE) pixelsHigh:(int)(YSIZE)
  461.         bitsPerSample:8 samplesPerPixel:3
  462.         hasAlpha:NO isPlanar:NO colorSpace:NX_RGBColorSpace
  463.         bytesPerRow:(int)(XSIZE*3) bitsPerPixel:0];
  464.     if (imageRep == nil) {
  465.         fprintf(stderr, "Couldn't create ImageRep.");
  466.         return NO;
  467.     }
  468.      
  469.     fd = open(name, O_CREAT | O_WRONLY | O_TRUNC, 0666);
  470.     if (!fd) { printf("Error opening save file %s.\n", name); }
  471.     stream = NXOpenFile(fd, NX_WRITEONLY);
  472.     if (stream == NULL) {
  473.         fprintf(stderr, "Error opening save stream.\n");
  474.         return NO;
  475.     }
  476.     [imageRep writeTIFF:stream usingCompression:NX_TIFF_COMPRESSION_LZW];
  477.     NXFlush(stream);
  478.     NXClose(stream);
  479.     close(fd);
  480.     return YES;
  481. }
  482.  
  483. - windowWillResize:sender toSize:(NXSize *)frameSize    // watch our window
  484. {
  485.     if (frameSize->width > (MAXXSIZE + 2)) frameSize->width = MAXXSIZE + 2;
  486.     if (frameSize->height > (MAXYSIZE + 32)) frameSize->height = MAXYSIZE + 32;
  487.     if (frameSize->width < 66) frameSize->width = 66;
  488.     if (frameSize->height < 96) frameSize->height = 96;
  489.     return self;
  490. }
  491.  
  492. - windowDidResize:sender                                // watch our window
  493. {
  494.     [saveMenuButton setEnabled:NO];
  495.     [reApplyButton setEnabled:NO];
  496.     return self;
  497. }
  498.  
  499.  
  500.  
  501. @end
  502.